Analyse Individual Dynos Captured with Stereo-seq¶
In [1]:
import stereo as st
import scanpy as sc
import pandas as pd
import numpy as np
import matplotlib as plt
import sys
import os
import matplotlib.pyplot as plt
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html from .autonotebook import tqdm as notebook_tqdm
Information on chip:
- horizontal dynos: Control + chir
- vertial dynos: DYMK + chir
Load in chip at bin20 then section different dynos for analysis¶
In [2]:
## Load into stereoexpdata object
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
spatial_key='spatial')
dyno_dm
Out[2]:
StereoExpData object with n_cells X n_genes = 13099 X 42359 bin_type: None offset_x = None offset_y = None cells: ['cell_name', 'total_counts', 'pct_counts_mt', 'n_genes_by_counts'] genes: ['gene_name'] cells_matrix = ['spatial'] Layers with keys: tl.result: []
In [3]:
ins = dyno_dm.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)
ins.show()
Out[3]:
In [4]:
dyno_dm_1 = dyno_dm.tl.filter_coordinates(min_x=9600,max_x=12600,min_y=12600,max_y=14600)
ins = dyno_dm_1.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)
dyno_dm_1.tl.filter_cells(
min_gene=200,
min_n_genes_by_counts=3,
pct_counts_mt=5,
inplace=True
)
dyno_dm_1 = st.io.stereo_to_anndata(dyno_dm_1,flavor='seurat',
output='dyno_dm_1.h5ad')
ins.show()
[2025-04-14 14:47:10][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates... [2025-04-14 14:47:10][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0054s. [2025-04-14 14:47:13][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells... [2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0404s. [2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident']. [2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] . [2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] . [2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info. [2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata. [2025-04-14 14:47:14][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_1.h5ad
Out[4]:
In [5]:
## Load into stereoexpdata object
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
spatial_key='spatial')
dyno_dm_2 = dyno_dm.tl.filter_coordinates(min_x=10000,max_x=11500,min_y=5700,max_y=9300)
ins = dyno_dm_2.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)
dyno_dm_2.tl.filter_cells(
min_gene=200,
min_n_genes_by_counts=3,
pct_counts_mt=5,
inplace=True
)
dyno_dm_2 = st.io.stereo_to_anndata(dyno_dm_2,flavor='seurat',
output='dyno_dm_2.h5ad')
ins.show()
[2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates... [2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0013s. [2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells... [2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0373s. [2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident']. [2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] . [2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] . [2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info. [2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata. [2025-04-14 14:47:15][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_2.h5ad
Out[5]:
In [6]:
## Load into stereoexpdata object
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
spatial_key='spatial')
dyno_dm_3 =dyno_dm.tl.filter_coordinates(min_x=13600,max_x=16600,min_y=16600,max_y=19600)
ins = dyno_dm_3.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)
dyno_dm_3.tl.filter_cells(
min_gene=200,
min_n_genes_by_counts=3,
pct_counts_mt=5,
inplace=True
)
dyno_dm_3 = st.io.stereo_to_anndata(dyno_dm_3,flavor='seurat',
output='dyno_dm_3.h5ad')
ins.show()
[2025-04-14 14:47:16][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates... [2025-04-14 14:47:16][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0014s. [2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells... [2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0761s. [2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident']. [2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] . [2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] . [2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info. [2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata. [2025-04-14 14:47:17][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_3.h5ad
Out[6]:
In [7]:
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
spatial_key='spatial')
dyno_dm_4 = dyno_dm.tl.filter_coordinates(min_x=9600,max_x=12600,min_y=16000,max_y=18600)
dyno_dm_4
dyno_dm_4.tl.filter_cells(
min_gene=200,
min_n_genes_by_counts=3,
pct_counts_mt=5,
inplace=True
)
ins = dyno_dm_4.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)
dyno_dm_4 = st.io.stereo_to_anndata(dyno_dm_4,flavor='seurat',
output='dyno_dm_4.h5ad')
ins.show()
[2025-04-14 14:47:18][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates... [2025-04-14 14:47:18][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0015s. [2025-04-14 14:47:18][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells... [2025-04-14 14:47:18][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0763s. [2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident']. [2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] . [2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] . [2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info. [2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata. [2025-04-14 14:47:19][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_4.h5ad
Out[7]:
In [8]:
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
spatial_key='spatial')
dyno_dm_5= dyno_dm.tl.filter_coordinates(min_x=11700,max_x=13200,min_y=3200,max_y=10150)
dyno_dm_5.tl.filter_cells(
min_gene=200,
min_n_genes_by_counts=3,
pct_counts_mt=5,
inplace=True
)
ins = dyno_dm_5.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)
dyno_dm_5 = st.io.stereo_to_anndata(dyno_dm_5,flavor='seurat',
output='dyno_dm_5.h5ad')
ins.show()
[2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates... [2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0013s. [2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells... [2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0363s. [2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident']. [2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] . [2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] . [2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info. [2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata. [2025-04-14 14:47:20][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_5.h5ad
Out[8]:
In [9]:
dyno_dm = st.io.read_ann_h5ad("/group/tran3/NGS/raw-data/BGI_STOmics_grant-08_2023/SAW_outputs/dyno_dymk_v2/051.cellcluster/A03390E1.cell.cluster.h5ad",
spatial_key='spatial')
dyno_dm_6= dyno_dm.tl.filter_coordinates(min_x=13000,max_x=15000,min_y=3400,max_y=10500)
ins = dyno_dm_6.plt.interact_spatial_scatter(width=500, height=500, poly_select=True)
dyno_dm_6.tl.filter_cells(
min_gene=200,
min_n_genes_by_counts=3,
pct_counts_mt=5,
inplace=True
)
dyno_dm_6 = st.io.stereo_to_anndata(dyno_dm_6,flavor='seurat',
output='dyno_dm_6.h5ad')
ins.show()
[2025-04-14 14:47:21][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_coordinates... [2025-04-14 14:47:21][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_coordinates end, consume time 0.0013s. [2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][st_pipeline][41][INFO]: start to run filter_cells... [2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][st_pipeline][44][INFO]: filter_cells end, consume time 0.0804s. [2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][940][INFO]: Adding sample in adata.obs['orig.ident']. [2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][943][INFO]: Adding data.position as adata.obsm['spatial'] . [2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][948][INFO]: Adding data.position as adata.obs['x'] and adata.obs['y'] . [2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][1127][INFO]: Rename QC info. [2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][1189][INFO]: Finished conversion to anndata. [2025-04-14 14:47:22][Stereo][3315164][MainThread][140157106832512][reader][1193][INFO]: Finished output to dyno_dm_6.h5ad
Out[9]:
In [10]:
import anndata as ad
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
from shapely.geometry import Polygon, Point
#import decoupler as dc
import scanpy as sc
from sklearn.neighbors import NearestNeighbors
from scipy import stats
import squidpy as sq
#from pydeseq2.dds import DeseqDataSet, DefaultInference
#from pydeseq2.ds import DeseqStats
In [11]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA
adata1=ad.read('dyno_dm_1.h5ad')
# Extract spatial coordinates
points = adata1.obsm['spatial']
# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)
# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
point = Point(row['x'], row['y'])
return point.distance(polygon.boundary)
# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
compute_distance_to_boundary, polygon=hull_polygon, axis=1
)
# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()
# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']
# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0] # The first principal component represents the longer axis
# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis
# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj
# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3
# Classification function
def classify_region(projection, threshold_1, threshold_2):
if projection <= threshold_1 or projection >= threshold_2:
return 'pole'
else:
return 'center'
# Apply classification on filtered points
adata_filtered.obs['region'] = [
classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')
# Visualize results
sq.pl.spatial_scatter(
adata_filtered,
size=10,
library_id="spatial",
figsize=(5, 5),
frameon=False,
shape=None,
color=["region"]
)
adata_filtered.write('dyno1.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored _cax = scatter(
In [8]:
from matplotlib.legend import Legend
In [12]:
adata_filtered=ad.read('dyno1.h5ad')
sq.pl.spatial_scatter(
adata_filtered,
size=10,
library_id="spatial",
figsize=(5, 5),
frameon=False,
shape=None,
color=["region"],title=''
)
fig = plt.gcf()
# Remove the legend
for ax in fig.get_axes():
for obj in ax.get_children():
if isinstance(obj, Legend):
obj.remove()
# **Manually Find and Remove the Colorbar**
for ax in fig.axes:
if "colorbar" in str(ax): # Check if this is the colorbar axis
fig.delaxes(ax)
# **Second Attempt: Remove Extra Axes (Colorbars are usually last axes added)**
while len(fig.axes) > 1: # Assuming 1 main plot, remove extras
fig.delaxes(fig.axes[-1]) # Delete last added axis (likely colorbar)
plt.savefig('/group/tran3/mnucera/STOmics_crossorganoid_project/STOmics_notebooks/figures_stomics/fig4A.svg',dpi=300)
plt.savefig('/group/tran3/mnucera/STOmics_crossorganoid_project/STOmics_notebooks/figures_stomics/fig4A.svg',dpi=300)
plt.show()
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored _cax = scatter(
In [13]:
adata_filtered=ad.read('dyno5.h5ad')
sq.pl.spatial_scatter(
adata_filtered,
size=10,
library_id="spatial",
figsize=(5, 5),
frameon=False,
shape=None,
color=["region"],title=''
)
fig = plt.gcf()
# Remove the legend
for ax in fig.get_axes():
for obj in ax.get_children():
if isinstance(obj, Legend):
obj.remove()
# **Manually Find and Remove the Colorbar**
for ax in fig.axes:
if "colorbar" in str(ax): # Check if this is the colorbar axis
fig.delaxes(ax)
# **Second Attempt: Remove Extra Axes (Colorbars are usually last axes added)**
while len(fig.axes) > 1: # Assuming 1 main plot, remove extras
fig.delaxes(fig.axes[-1]) # Delete last added axis (likely colorbar)
plt.savefig('/group/tran3/mnucera/STOmics_crossorganoid_project/STOmics_notebooks/figures_stomics/fig4B.svg',dpi=300)
plt.show()
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored _cax = scatter(
In [3]:
def create_custom_legend(labels, colors, filename="legend.svg", dpi=300):
fig, ax = plt.subplots()
ax.axis("off") # Hide axes
# Create legend handles
handles = [plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=color, markersize=10)
for color in colors]
# Add legend
legend = ax.legend(handles, labels, loc="center", frameon=False)
# Save figure with tight bounding box to remove excess whitespace
fig.savefig(filename, dpi=dpi, bbox_inches="tight", transparent=True)
plt.show()
# Example usage
labels = ["Center","Poles"]
colors = ['blue','orange']
create_custom_legend(labels, colors, filename="/group/tran3/mnucera/STOmics_crossorganoid_project/STOmics_notebooks/figures_stomics/dyno_legend.svg")
In [16]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA
adata1=ad.read('dyno_dm_2.h5ad')
# Extract spatial coordinates
points = adata1.obsm['spatial']
# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)
# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
point = Point(row['x'], row['y'])
return point.distance(polygon.boundary)
# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
compute_distance_to_boundary, polygon=hull_polygon, axis=1
)
# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()
# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']
# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0] # The first principal component represents the longer axis
# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis
# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj
# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3
# Classification function
def classify_region(projection, threshold_1, threshold_2):
if projection <= threshold_1 or projection >= threshold_2:
return 'pole'
else:
return 'center'
# Apply classification on filtered points
adata_filtered.obs['region'] = [
classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')
# Visualize results
sq.pl.spatial_scatter(
adata_filtered,
size=10,
library_id="spatial",
figsize=(5, 5),
frameon=False,
shape=None,
color=["region"]
)
adata_filtered.write('dyno2.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored _cax = scatter(
In [17]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA
adata1=ad.read('dyno_dm_3.h5ad')
# Extract spatial coordinates
points = adata1.obsm['spatial']
# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)
# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
point = Point(row['x'], row['y'])
return point.distance(polygon.boundary)
# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
compute_distance_to_boundary, polygon=hull_polygon, axis=1
)
# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()
# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']
# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0] # The first principal component represents the longer axis
# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis
# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj
# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3
# Classification function
def classify_region(projection, threshold_1, threshold_2):
if projection <= threshold_1 or projection >= threshold_2:
return 'pole'
else:
return 'center'
# Apply classification on filtered points
adata_filtered.obs['region'] = [
classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')
# Visualize results
sq.pl.spatial_scatter(
adata_filtered,
size=10,
library_id="spatial",
figsize=(5, 5),
frameon=False,
shape=None,
color=["region"]
)
adata_filtered.write('dyno3.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored _cax = scatter(
In [18]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA
adata1=ad.read('dyno_dm_4.h5ad')
# Extract spatial coordinates
points = adata1.obsm['spatial']
# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)
# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
point = Point(row['x'], row['y'])
return point.distance(polygon.boundary)
# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
compute_distance_to_boundary, polygon=hull_polygon, axis=1
)
# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()
# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']
# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0] # The first principal component represents the longer axis
# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis
# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj
# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3
# Classification function
def classify_region(projection, threshold_1, threshold_2):
if projection <= threshold_1 or projection >= threshold_2:
return 'pole'
else:
return 'center'
# Apply classification on filtered points
adata_filtered.obs['region'] = [
classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')
# Visualize results
sq.pl.spatial_scatter(
adata_filtered,
size=10,
library_id="spatial",
figsize=(5, 5),
frameon=False,
shape=None,
color=["region"]
)
adata_filtered.write('dyno4.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored _cax = scatter(
In [19]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA
adata1=ad.read('dyno_dm_5.h5ad')
# Extract spatial coordinates
points = adata1.obsm['spatial']
# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)
# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
point = Point(row['x'], row['y'])
return point.distance(polygon.boundary)
# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
compute_distance_to_boundary, polygon=hull_polygon, axis=1
)
# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()
# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']
# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0] # The first principal component represents the longer axis
# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis
# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj
# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3
# Classification function
def classify_region(projection, threshold_1, threshold_2):
if projection <= threshold_1 or projection >= threshold_2:
return 'pole'
else:
return 'center'
# Apply classification on filtered points
adata_filtered.obs['region'] = [
classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')
# Visualize results
sq.pl.spatial_scatter(
adata_filtered,
size=10,
library_id="spatial",
figsize=(5, 5),
frameon=False,
shape=None,
color=["region"]
)
adata_filtered.write('dyno5.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored _cax = scatter(
In [ ]:
In [20]:
import numpy as np
import scanpy as sc
import squidpy as sq
from scipy.spatial import ConvexHull
from shapely.geometry import Point, Polygon
from sklearn.decomposition import PCA
adata1=ad.read('dyno_dm_6.h5ad')
# Extract spatial coordinates
points = adata1.obsm['spatial']
# Compute the convex hull
hull = ConvexHull(points)
hull_points = points[hull.vertices]
hull_polygon = Polygon(hull_points)
# Function to compute distance from point to polygon boundary
def compute_distance_to_boundary(row, polygon):
point = Point(row['x'], row['y'])
return point.distance(polygon.boundary)
# Apply function to all spots
adata1.obs['distance_to_boundary'] = adata1.obs.apply(
compute_distance_to_boundary, polygon=hull_polygon, axis=1
)
# Define threshold to remove perimeter points
edge_threshold = adata1.obs['distance_to_boundary'].quantile(0.15)
adata_filtered = adata1[adata1.obs['distance_to_boundary'] > edge_threshold].copy()
# Extract filtered spatial coordinates
filtered_points = adata_filtered.obsm['spatial']
# Perform PCA to find the principal axis (longer axis) on filtered points
pca = PCA(n_components=2)
pca.fit(filtered_points)
principal_axis = pca.components_[0] # The first principal component represents the longer axis
# Project filtered points onto the principal axis
projected_distances = filtered_points @ principal_axis
# Find min and max projections (endpoints of the longer axis)
min_proj, max_proj = projected_distances.min(), projected_distances.max()
total_length = max_proj - min_proj
# Define threshold to split into three sections
threshold_1 = min_proj + total_length / 3
threshold_2 = min_proj + 2 * total_length / 3
# Classification function
def classify_region(projection, threshold_1, threshold_2):
if projection <= threshold_1 or projection >= threshold_2:
return 'pole'
else:
return 'center'
# Apply classification on filtered points
adata_filtered.obs['region'] = [
classify_region(proj, threshold_1, threshold_2) for proj in projected_distances
]
adata_filtered.obs['region'] = adata_filtered.obs['region'].astype('category')
# Visualize results
sq.pl.spatial_scatter(
adata_filtered,
size=10,
library_id="spatial",
figsize=(5, 5),
frameon=False,
shape=None,
color=["region"]
)
adata_filtered.write('dyno6.h5ad')
/group/neur6/singularity_imgs/.conda/envs/st2025/lib/python3.8/site-packages/squidpy/pl/_spatial_utils.py:955: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored _cax = scatter(
In [ ]: